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*^ Abstract Direct Monte Carlo simulations of the Enskog- rived in cases (a) and (b) [T. P. C. van Noije and M. H. 
D 

g Boltzmann equation for a spatially uniform system of smoothErnst, Gran. Matt. 1, 57 (1998)] shows a good agreement. 

+^ ■ inelastic spheres are performed. In order to reach a steady In contrast to these two cases, the deviation from the 

> \ 

^ . state, the particles are assumed to be under the action of Maxwell-Boltzmann distribution is not well represented 
+-* . 

g an external driving force which does work to compensate by Sonine polynomials in case (c), even for low dissipa- 
for the collisional loss of energy. Three different types of tion. In addition, the high energy tail exhibits an under- 

S3 : 

Q external driving are considered: (a) a stochastic force, (b) population effect in this case. 

O ' 

■ a deterministic force proportional to the particle velocity 



. and (c) a deterministic force parallel to the particle ve- 

m ; 

C^l ' locity but constant in magnitude. The Enskog-Boltzmann 

m : 

^| equation in case (b) is fully equivalent to that of the homo- 



o 
o 



geneous cooling state (where the thermal velocity mono- 
■ tonically decreases with time) when expressed in terms of 



1 

Introduction 

Most of the recent studies of rapid granular flow Q are 
based on the Enskog equation for the velocity distribution 



the particle velocity relative to the thermal velocity. Com- 
. function /(r, v, t) of an assembly of inelastic hard spheres 

parison of the simulation results for the fourth cumulant 



o 
o 



X 



and the high energy tail with theoretical predictions de- 



ll. In the special case of a spatially uniform state, the 
Enskog equation reads 
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function at contact and / is the collision operator. In 
Eq. (||), d is the dimensionality of the system, a is the 
diameter of the spheres, V12 = vi — V2 is the relative 
velocity of the colliding particles, <r is a unit vector di- 
rected along the line of centers from the sphere 1 to the 
sphere 2, is the Heaviside step function and a < 1 
is the coefficient of normal restitution, here assumed to 
be constant. In addition, (v",v 2 ') are the precollisional 
velocities yielding (vi, V2) as the postcollisional ones, i.e. 
v i 2 = v i,2T ^(l + a _1 )(vi2 -<x)<x. Except for the presence 
of the factor x, which accounts for the increase of the colli- 
sion frequency due to excluded volume effects, the Enskog 
equation for uniform states, Eq. ([l]), becomes identical 
with the Boltzmann equation. 

In the case of elastic particles (a = 1) and in the ab- 
sence of external forcing (T = 0), it is well known that the 
long-time solution of Eq. ([!]) is the Maxwell-Boltzmann 
equilibrium distribution function, /(v,i) — > nVg d (j)(v/vQ), 
0(c) = ii~ d / 2 e~ c , where n is the number density, vq = 
(2(w 2 )/(i) 1 / 2 is the thermal velocity and c = v/i?o is the 
reduced velocity. On the other hand, if the particles are 
inelastic (a < 1) and T = 0, a steady state is not possi- 
ble in uniform situations since, due to the dissipation of 
energy through collisions, the thermal velocity vo(t) de- 
creases monotonically with time. Regardless of the initial 
uniform state, the solution of Eq. ([[]) tends to the so-called 
homogeneous cooling state characterized by the 

fact that the time dependence occurs only through the 
thermal velocity vo(t): /(v, i) — > nv^ {t)f iy / Voiff) . In 
addition, /(c) deviates from a Maxwellian, /(c) 7^ 0(c), 



as measured by the fourth cumulant 
where 

( C P) = J dccPf(c). (4) 

By expanding /(c) /0(c) in a set of Sonine polynomials 
{S p (c 2 )} and neglecting the terms beyond p = 2, van Noije 
and Ernst have estimated the value of a^. 

( s 16(1 -«)(!- 2a 2 ) 

a2(a> ~ 9 + 24d - a(41 - 8d) + 30(1 - a)a 2 ' 1 ' 

The above expression corrects an algebraic error in a pre- 
vious calculation of a-i in the three-dimensional case 
According to Eq. (||), 02 changes sign at a = l/s/2 ~ 0.71. 
By using the same method, Garzo and Dufty have re- 
cently extended the evaluation of 02 to a binary mixture 
of hard spheres. The accuracy of Eq. (g) has been quan- 
titatively confirmed by Brey et al. [0] from Monte Carlo 
simulations of the Boltzmann equation for hard spheres 
(d = 3) in the range 0.7 < a < 1. As a complementary 
measure of the departure of /(c) from 0(c), Esipov and 
Poschel II and van Noije and Ernst || have analyzed 
the high energy tail of the distribution function and have 
found an asymptotic behavior of the form 

log/(c)~-c, (6) 

in contrast to log 0(c) ~ — c 2 . The high energy tail (^) has 
been confirmed by simulations in the case of hard disks 
(d = 2) §. 

In order to reach a steady state, energy injection is 
needed to compensate for the energy dissipated through 
collisions. This can be achieved by vibration of vessels 



JTo[ or in fluidized beds [jO]. The same effect can be ob- where £ is a positive constant. In this case, 

tained by means of external driving forces acting locally rt , <. > d . . 

66 J7(vi) = C— • [vi/(vi)j . (13) 

on each particle ]l2] ] . Borrowing a terminology frequently 
used in nonequilibrium molecular dynamics of elastic par- 
ticles jUI, we will call this type of external forces "ther- 
mostats" . In general, the equation of motion for a particle 
i is then 



<9vi 

It is interesting pointing out that the Enskog-Boltzmann 
equation (|l|) for the above Gaussian thermostat force is 
formally identical with the equation for the homogeneous 
cooling state (i.e. with T = 0) when both equations are 
expressed in terms of the reduced distribution /(c) (see 
mVi = F « + F * ' ( 7 ) Sect. |). As a consequence, the results (§) and (|) apply 

where m is the mass of a particle, F™ 11 is the force due to to this thermostatted case as well. 

collisions and F- h is the thermostat force. Williams and The differences between Eqs. @ and © and between 
MacKintosh (12) introduced a stochastic force assumed to ® and (||) illustrate the influence of the thermostat force 
have the form of a Gaussian white noise: on tlie departure of the steady-state distribution function 

from the Maxwell-Boltzmann distribution. In the case of 



<Ff(i)> = 0, (F?(t)Ff(t')) = \m 2 e 5 ij S(t-t'), (8) 

where I is the d x d unit matrix and £q represents the 
strength of the correlation. The corresponding operator 
T appearing in Eq. ([|) has a Fokker-Planck form |j| : 



the stochastic force, Eq. (^|), /(c) is closer to <j>(c) than in 



the case of the Gaussian force, Eq. (13), since in the former 
case ci2 and the high energy overpopulation are smaller 
than in the latter case. Of course, other types of ther- 



£ 2 f 9 V 

■7\f( v i) = — ~2 ( ~q^~ J /( v i)- (9) mostats are also possible. For example, a different choice 

fr it • • , t-, , rci , ij-jii. i j • , for a deterministic thermostat is 

Van JNoije and Ernst [bfl have studied the stationary solu- 
tion of the uniform equation ([!]) with the thermostat (^). F* h = mgVi, (14) 

They have found for the coefficient a-i defined by Eq. (0) i — i „„ -1^1^ e t-> ,177^ • 

^ where v, = w l /v t . While the Gaussian force, Eq. (|12|), is 

proportional to the velocity of the particle, Eq. (R4h corre- 
cts (a) ~ — — ^- ?^_) 1 iq\ sponds to a force that is parallel to the direction of motion 

21 ' ~ 73 + 56d - 3a(35 + 8<i) + 30(1 - a)a 2 K ' 

The high energy tail is || 

log/(c)~-c 3 / 2 . (11) 

Of course, deterministic thermostats can also be used. 
For instance, the use of Gauss's principle of least con- 
straint leads to the thermostat force |D| 



but constant in magnitude. The corresponding operator T 



is 



mCvi, (12) 



•77(vi) = ^ { ^ ■ [vi/( Vl )] - /( Vl )| • (15) 

The aim of this paper is to present direct Monte Carlo 
simulations of Eq. (Q) with the three choices for the ther- 
mostat, Eqs. (|), (|f) and ( fl5| ) . In the cases of the stochas- 
tic and the Gaussian thermostats, we will confirm the tails 



(|Tl|) and (H) and will check the accuracy of the estimates where 

(llCl) and (H). In the latter case, however, we will see that ~ f f . A ^ 

^ I[ Cl \f,f}=j dc 2 J dff©(ci 2 -ff)(cia-ff) 

a better agreement with simulation results for a < 0.5 is r „ „ ~ ~ n 

x [a- 2 /(c?)/(40 - /(ci)/(c 2 )j . (17) 

obtained if an estimate slightly different from (|5|) is used. 

The simulation results corresponding to the non-Gaussian The reduced operator T for the stochastic [Eq. ®] , Gaus- 

thermostat @ show that, in contrast to what happens in sian [ Eq - ©1 and non-Gaussian [Eq. @] thermostats, 



d _ x g/(ci) 
1 9 Cl 



(18) 



the two previous cases, ai remains negative for all a.. This 1S 

feature is qualitatively captured by an estimate derived ~~» . £n -(d-i) 9 

2v£ncr d ~ i oci 

from a Sonine approximation. In this problem, however, 

the Soninc polynomials do not constitute a good set for the . _ 

Ff(c x ) = — Vrcr (d - 1J k/( Cl )j , (19) 
expansion of f(c)/<f>(c) and, consequently, the estimate is vgna oc\ i J 



-(d-i) 3 



not quantitatively good. Besides, the high energy tail is of 

.T/(ci)= „ g— cT^-^-Jl- cf-VCci) , (20) 
the form log /(c) ~ — c , but with a coefficient different v^na 11 1 acj L J 

from that of the Maxwell-Boltzmann distribution. respectively. In Eqs. ©-© we have already taken into 

The organization of this paper is as follows. The theo- account that the distribution function must be isotropic 
retical analysis is reviewed in Sect. |. The computer simu- in the steady state. Equation © with the term @ is 
lation method employed to solve numerically the uniform fully equivalent to Eq. (10) of Ref. g, the latter being 
Enskog-Boltzmann equation is described in Sect. § The derived in the context of the homogeneous cooling state, 
results are presented and compared with the theoretical This formal equivalence between the free evolving state 
predictions in Sect. |. The paper ends with a summary and the one controlled by a Gaussian external force is 
and discussion in Sect. ^. also present in the case of elastic particles interacting via 

arbitrary power-law potentials in homogeneous situations 
|L4[ or via the Maxwell potential in the uniform shear flow 

2 S 
Theoretical predictions 

2.1 

In the steady state, the Enskog-Boltzmann equation ([!]) Stochastic thermostat 

can be expressed in terms of the reduced velocity distri- „ . , , . , . . P 

l*or the sake ot completeness, we summarize now some ot 

bution function /(c) as results obtained in Ref. In order to characterize the 

deviation of /(c) from 0(c) by means of the cumulant (||), 

^/(ci) = x7[ Cl |/,7], (16) 

it is uselul to consider the hierarchy ot moment equations. 



Multiplying both sides of Eq. (|16|) by of and integrating 
over Ci , we get 



,(!) - 



(0) 



Mp = M2- 



p(p + d - 2) 
2d 



„P~2\ 



(21) 



for the stochastic thermostat, where we have defined 



dcc*7[c|/,/]. 



(22) 



In Eq. (|21j) we have taken into account the normalization 
condition (c°) = 1, so that /12 = ^Co/' i; oX ncrd_1 - I n the 
special case of p — 4, Eq. ( pi] ) becomes 

m = (d+2)fi 2 , (23) 

where we have used the fact that, by definition, (c 2 ) = 
d/2. Equations ( p2| ) and (^3|) are still exact. To get an 
approximate expression for 02(0;), three steps will be taken 
H. First, we assume that / can be well described by the 
simplest Sonine approximation, at least for the velocities 
relevant to the evaluation of 02- Thus, 

/(c) ~ 0(c) [l + a2 S 2 (c 2 )] , (24) 

where 



The approximation (J2J) is justified by the fact that 02 is 
expected to be small. The second step consists of inserting 
Eq. (^4|) into Eq. (^2|) and neglecting terms nonlinear in 
Q2- For /Z2 and /Z4 the results are B 



with 



t<2 



(0) _ " f-, „.2\ 



(26) 



v / 2r(d/2) 



(l-« 2 ), 



(27) 



Mi 0) 



^2 = y^ 2 



— (l0ci + 39+ 10c 



d- 1 



1 - a 



(0) 

/4 ■ 



(28) 
(29) 

(30) 



In the third step, the approximations ( pq ) with p = 2 and 4 
are inserted into the exact equation ( p3| ) and 02 is obtained 
from the resulting linear equation: 

^-(d + 2)^ 



a-2 



nP-(d+2)£ y 



(31) 



This is the result derived by van Noije and Ernst ||, Eq. 
( |l0| ) . It must be pointed out that a certain degree of am- 
biguity is present in this last step. For instance, if Eq. (|2J 
were written as /X4//Z2 = d + 2, we could expand the ratio 

M4/V2 in powers of a 2 and neglect nonlinear terms to find 
^-(d + 2)^ 



('2 



4(1 - a)(l - 2a 2 ) 



19+ 14d - 3a(9 + 2<i) + 6(1 -a)a 2 ' ^ 
However, since a 2 is indeed small (\a 2 \ < 0.1), Eqs. © 

and ( [32| ) give practically identical results, the maximum 

deviation being less than about 0.001. 

Now we consider the high energy tail. In general, the 

collision integral can be decomposed into a gain and a loss 

term: 7[ci|/, /] = ? 9 [ci|/, J] - Tt[cx\f, /]. For large c x the 

loss term can be approximated as 

%[c 1 \fj]=p 1 J cic 2 012/(00/(02) 

« PicJ(ci), (33) 

where (3\ is defined by Eq. (|S4|). Let us assume that for 
large velocities the gain term is negligible versus the loss 
term, i.e. 



hm 



/i[ci|/,/] 



= 0. 



(34) 



In that case, the Enskog-Boltzmann equation for the stochas- 
tic thermostat becomes 



2d dc 



,-i5/(c) 



dc 



Ac/(c). (35) 



The solution of this equation for large c is 



/(c) w Aexp (-Ac 3/2 



A 



2 /2d/3i 



1/2 



(36) 



3 V M2 

where If is an undetermined constant. By arguments given 
in Ref. ||, it can be seen that the result ( |36] ) is indeed 
consistent with the assumption (|34|) . Equation (|3^) shows 
an overpopulation with respect to the Maxwell-Boltzmann 
tail. On the other hand, as a — > 1, the amplitude A di- 
verges as (1 — a) -1 / 2 , thus indicating that the overpopu- 
lation effect is restricted to larger and larger energies in 
the limit a — ► 1. 

2.2 

Gaussian thermostat 

In the case of the deterministic Gaussian thermostat, Eq. 
(p^), the moment equation is 



where now /i 2 = d£/voxno~ . If we set p = 4 



(37) 



(j,i = (d+ 2)(1 + a 2 )M2> 



(38) 



connection with the approximation ( 26 ) . If we rewrite (|38 
as n±/ H2 = (d + 2)(1 + 02) and neglect nonlinear terms, 
the resulting a 2 is fairly close to Eq. (|39|). On the other 
hand, if we start from [14/(1 + a 2 ) = (d+ 2)/i 2 , the result 



is 



02 



M4 



(d+2) M2 ( 



16(l-a)(l-2a 2 ) 



25 + 24rf-a(57-8d)-2(l-a)a 2 ' < - 4 °" ) 
The estimates @ and (|^) practically coincide in the re- 
gion 0.5 < a < 1. However, they visibly separate for larger 
dissipation. In the interval < a < 0.3, the values given 
by Eq. (§) are 12%-2Q% (d = 3) or 18%-28% (d = 2) 
larger than those given by Eq. (fiof). As we will see later, 
the simulation results indicate that Eq. (]!(]) is a better 
estimate than Eq. (|J). 

For large c the Enskog-Boltzmann equation becomes 



tl c -(d-i)l 
d dc 



c d f(c) 



-/3ic/(c), 



(41) 



where we have used Eqs. ( |33| ) and (|34|). Its solution is 

/(c) « Aexp (-Ac) , A = — . (42) 

M2 

Again, this result is seen to be consistent with ( |34j ) [9. 
Equation ([|2]) indicates an overpopulation effect even larger 
than with the stochastic thermostat. 



2.3 



where we have made use of Eq. (||). Substituting the ap- Non-Gaussian thermostat 

proximation © and neglecting terms nonlinear in a 2 , we Now we consider the deterministic non-Gaussian thermo- 



get 



a 2 



^-(d + 2)^ 



(39) 



^-(d+2)(^+^y 

which is the same as Eq. (||). There exists again some 
arbitrariness about the use of the exact equation ([58]) in 



stat (|14|), represented by the operator (J20[) . To the best of 
our knowledge, this external force has not been analyzed 
before. The corresponding moment equation is 



P (c 



(43) 



where fi 2 = 2g(c) jv^xna . In particular, case, fj, 2 and /X4 would be given by Eqs. (|26|)-(|30D and 



»^C + W/ 1 + « |. (49) 



In contrast to the two previous cases, now the even colli- 

sional moments \i v are coupled to the odd moments (c^" 1 ), ^(d/2) 
and vice versa. In terms of the energy variable e = c 2 , Inserting this into Eq. Q) and neglecting nonlinear terms 
this means that the integer collisional moments are cou- we get 

pled to the half-integers energy moments. This is related — (d + 1)^ 2 °^ 

a 2 ~ 



- (d+l)(u {1) + a (0) /2) 

to the fact that the force (1141) is singular at e = 0. As ^ K A 2 2 'J 

16(l-a)(l + 2a 2 ) 

a consequence, while /(c) is expected to be close to the 63 + 40d — a(95 + 8d) + 30(1 — a)a 2 

Maxwellian 0(c), the ratio /(c)/</>(c) is singular at e = While in the cases of the stochastic thermostat, Eqs. @ 
and thus it is not well represented by an expansion in or ©, and the Gaussian thermostat, Eqs. (§) or ©, 
{5 p (e)}. To be more precise, let us define the function the cumulant a 2 changes from negative to positive values 



A(c) by the equation at a ~ 0.71, Eq. ( |5C| ) indicates that 02 remains negative 

_ in the case of the non-Gaussian thermostat. We will see in 

7(c) = 0(c) [l + a 2 4(c)]. (45) 

Sect. that our computer simulations confirm this feature. 
Therefore, At a quantitative level, however, the estimate ( |50| ) is about 



dc<f){c)A{c) = / dc0(c)cM(c) = 0, (46) 



20% too small in magnitude. 

To analyze the high energy tail, let us assume for the 



, , . . a ., . did + 2) , . moment the validity of (03) , so that the Enskog-Boltzmann 

dc (f>(c)c 4 A(c) = — -. (47) V- 2 " 



4 

The polynomial ^(c 2 ) verifies the above equalities. As a 



c -( d -i) d 

matter of fact, A(c) ~ ^(c 2 ) in the cases of the ther- 2(c) dc 



equation can be replaced by 

c^fic) a-Ac/(c), (.VI) 



mostats (|18|) and (|19|). This is not so, however, in the case wnose solution for large c is 
of (p0|), even in the limit of low dissipation. As we will see 
in Sect. |, dA(c)/dc\ c=0 ^ 0, what indicates that A(c) 



/(c) «^exp(-A'c 2 ), A' = ^-. (52) 



is essentially different from a polynomial in c 2 . All of this According to (52), /(c) has a Maxwellian tail that is un- 

complicates the evaluation of a 2 . Nevertheless, since A(c) derpopulated with respect to the Maxwell-Boltzmann dis- 

and S 2 (c 2 ) share the moments of degrees 0, 2 and 4 [cf. tribution 0(c), since the amplitude A' ~ \/2/(l — a 2 ) is 

Eqs. ( |46| ) and (f47j)] , we can expect to obtain a crude esti- larger than 1. But now we get an unphysical result: the un- 

mate of a 2 by assuming that in the calculation of (c), (c 3 ), derpopulation effect increases as one approaches the elas- 

\i 2 and /Z4 we can replace A(c) by S 2 (c 2 ). If that were the tic limit, since A' — > 00 as a — > 1. The solution to this 



paradox lies in the fact that the assumption (|3j) is not where vq(0) is an arbitrary initial thermal velocity. To 
justified in this case. Let us assume instead that the gain enforce a vanishing initial total momentum, the velocity 
and loss term are comparable, namely of every particle is subsequently subtracted by the amount 

/ 9 [ci|/J]«7/W(ci), (53) JV ~ 1 ^ V *(°)' 

The velocities are updated from time t to time t + h, 

where 7 < 1 is an unknown function of a. According to 

where the time step h is much smaller than the mean 

this, Eq. (|5^ ) is replaced by 

free time, by following two successive stages: collisions 
/(c) ps K cxp (—Ac ) , A = A (1 — 7). (54) anc j f ree streaming. In the collision stage, a sample of 

On physical grounds we expect that A -> 1 when a -f 1, l JVw ™' 1 P airs is chosen at random with cquiprobability, 

which implies that 7 -> 1- ^2(1 - a) in that limit. As will where ^ is an upper bound estimate of the probability 

be shown in Sect. |, comparison with simulation results that a particle collides per unit of time. For each pair ij 

confirms a behavior of the form (Eil) belonging to this sample, the following steps are taken: (1) 



a given direction is chosen at random with equiprob- 

3 ability; (2) the collision between particles i and j is ac- 

Direct Simulation Monte Carlo method cepted with a probability equal to 0(vy • <Xtf)< 4 W a; niaxi 

The Direct Simulation Monte Carlo (DSMC) method de- where U V = (Wxn)| v<* ■ | ; if the collision is accepted, 

vised by Bird @ has proven to be a very efficient tool postcollisional velocities are assigned to both particles: 
to solve numerically the Boltzmann equation. The DSMC T H 1 + a X v y ' In the case that in 

method has been recently extended to the Enskog equa- one of the col li sions "ij > w m», the estimate of w max is 

tion and its application to inelastic particles is straight- updated as w max = . 

forward §,[[§]. Here we briefly describe the specific method In the free streaming stage the velocity of every par- 

we have used to solve the uniform Enskog-Boltzmann equa- tide is changed according to the thermostat force under 

tion (0) in the case of a three-dimensional system (d — 3). consideration: 

The velocity distribution function is represented by the 
velocities {v.;} of N "simulated" particles: 

i N 

/(v,t)-nl£«(v<(t)-v). (55) where 

i=l 

1 rt+h 

At the initial state the particles are assigned velocities Wj = / dt'F th (t'). (58) 

m J t 

drawn from a Maxwell-Boltzmann probability distribu- 

In the case of the stochastic thermostat, Eq. (IJ), one has 



(57) 



tion: 



n- 1 /(v,0)= 7 r- 3 / 2 V 3 (0)e- , ' 2/t '« 2(0) , (56) ( Wi ) = 0, (w iWj ) = \$h6 i:j . (59) 



Consequently, each vector Wj is randomly drawn from the In the above equations, V12 = vi — v 2 , V12 = j(vi 4- 
Gaussian probability distribution V2). Starting from the exact expression ( |o4| ) and using 



P(w) = (2vt;lhf 2 e- w2 > 2e ° h . (60) 



(p5|), we arrive at the following formula for the numerical 
computation of \i v : 



l'i> = v o P 1 ^p( v «' v i)- ( 67 ) 



1 j 



In the case of deterministic external forces the velocity 

increment Wj is assigned in a more direct way. If the ther- ' '' 11 N 

mostat is the Gaussian one, Eq. (112), ™ . . ,. , , • , 

1 Ihe prime m trie summation means that we restrict our- 

w . _ /gCft _ ^ v . (qi^ selves to N' pairs ij randomly chosen out of the total 

number N(N — l)/2 of pairs in the system. This allows us 

In the case of the non-Gaussian thermostat defined by Eq. 

to compute (c p ) and fi p with similar accuracy within rea- 

©» 

sonable computer times. Once the steady state is reached, 

Wj = gh (yi — k) , (62) the relevant quantities are subsequently averaged over M 

^ _ independent instantaneous values, 

where the vector k = N~~ Vj is introduced to preserve 

In our simulations we have typically taken N = 2 x 

the detailed conservation of momentum, i.e. w « = 0. 

10 5 , N' = 10 7 and M = 10 3 . Since the thermal velocity 

The moments of the distribution are simply obtained 

is not constant in the transient regime, we have taken a 

as 

time-dependent time step h — O.OIA/Wi), where A = 

1 ST- 

^ = nY, v *> ^ = ( yP )/ v ^ ( 63 ) (VTZvrxna 2 )- 1 is the mean free path. 

i—1 

1 /2 

where v = (2(v 2 )/3) ' . The evaluat ion of the collisional 
moments n P , p — 2 and 4, is more complicated. In the 4 
Appendix it is shown that Results 

-2 -0-1 f 1 f 1 n \ p> 1 \ By using the numerical method described in the previous 

M P = n 2 v Q p W dvi / dv 2 /(vi)/(v2)*p(vi,va), y 5 P 

section, we have computed the steady-state values of the 

(64) 

first few moments (c p ) and [i v . We have also evaluated 



where 



Tr(l-a 2 ) 



the reduced velocity distribution function /(c). As a test 

*a(vi,v 2 )= J vf 2) (65) 

of the accuracy of the simulations and also to check that 

the steady state has been reached, we compare in Table 

[j] the values of /14 obtained directly from Eq. (^7|) with 



^ , s n I 5(1 - a ) 2 ta2 

^ 4 (vi,V 2 ) = -V U { jj -«I2^12 



(l-a 2 )(2 + a 2 ) 4 



12 those given by Eqs. (E3j), (pq) or (|4J). The values corre- 



+(3-a)(l + a) 



(V12 ' V 12 ) 2 - -vf 2 ^12 



sponding to a Maxwell-Boltzmann distribution, \ are 
(66) also included in the table. We can observe that the direct 
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Table 1. Comparison of the values of /J4 obtained directly with those obtained indirectly from fi2- 



Stochastic Gaussian non-Gaussian 



a 




Hi 


Eq. @ 


/'4 




/'4 


Eq. © 


0.2 


10.925 


12.157 


12.155 


13.881 


13.881 


8.744 


8.750 


0.4 


9.812 


10.602 


10.600 


11.494 


11.488 


7.631 


7.631 


0.6 


7.797 


8.036 


8.038 


8.213 


8.217 


5.811 


5.810 


0.8 


4.638 


4.499 


4.503 


4.414 


4.412 


3.335 


3.333 



and indirect routes to the computation of /14 disagree less 
than 0.1% in all the cases. The difference between /14 and 
jif^ is a measure of the departure of /(c) from 4>{c). 

Now we present the results separately for each one of 
the three thermostats considered. 

4.1 

Stochastic thermostat 

The basic quantity measuring the deviation of the dis- 
tribution function from the Maxwell-Boltzmann distribu- 
tion is the cumulant a2, Eq. (||). Figure [j] shows the a- 
dependence of the simulation values of a 2 , (p.2 — a4°^ ) / > 
(/i4 — M4°'')// i 4 1 ^ an d the theoretical estimate (|l0|), first de- 
rived in Ref. ||. As said in Sect. |], the estimate ([52]) 
gives practically the same results as (|l0|) and therefore 
it is not plotted. The agreement between the simulation 
data and the theoretical prediction is excellent, thus in- 
dicating that the approximation ( p6| ) was justified. The 
above agreement indicates that the distribution function 
/(c) for thermal velocities is well represented by Eq. (pi|). 
To confirm this, the function A(c) defined by Eq. (|4^) is 
plotted in Fig. ^ for a = 0.5. The simulation curve agrees 
very well with the Sonine polynomial ^(c 2 ). It is worth 




-0.02 I ' 1 ' 1 ' 1 ' 1 ' 1 

0.0 0.2 0.4 0.6 0.8 1.0 

a 

Fig. 1. Plot of the simulation values of a2 (0)> (^2 — /-4° )/a4 
(A) and (in—y^f^ 1 )/ fj^ (v) versus a in the case of the stochas- 
tic thermostat. The solid line is the theoretical estimate, Eq. 




noting that this deviation from the Maxwell-Boltzmann 
distribution in the case of the stochastic thermostat could 
not be observed in recent two-dimensional molecular dy- 
namics simulations [ fl9| because the statistical accuracy 
was not high enough. 

The theoretical prediction for the asymptotic high en- 
ergy tail, Eq. (|36|), is much harder to confirm in the simu- 
lations since it involves a very small fraction of particles. 
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A(c) 




G(c) 
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Fig. 2. Plot of the simulation values of the function A(c) dc- pig 3 plot of the simulation values of the functi on G(c) de- 
fined by Eq. (^) for a = 0.5 in the case of the stochastic 
thermostat. The dashed line is the Sonine polynomial 5*2(c 2 ). 



Equation (|3q ) implies that 

lim G(c) = K = const, 



(68) 



where 



G(c) = e Ac " /(c). 



(69) 



The function G(c) is plotted (in logarithmic scale) in Fig. 
^ for a = 0.4 and a = 0.5. In both cases the values of A 
have been obtained from (|36| ) by using the simulation val- 
ues of fi 2 , which yields A ~ 1.99 (a = 0.4) and A ~ 2.11 
(a = 0.5). The figure is convincingly consistent with Eq. 
(||), where K ~ 1.3 and K ~ 2.2 for a = 0.4 and a = 0.5, 
respectively. Figure |^ also shows the corresponding func- 
tions G(c) obtained from Eq. ( |69| ) by replacing /(c) by 
the Maxwell-Boltzmann distribution 0(c). The overpopu- 
lation phenomenon for c > 2 is quite apparent. At c = 4, 
for instance, f/<f> ~ 8 for a = 0.4 and f /4> — 5 for a = 0.5. 



fined by Eq. ( J69D for a = 0.4 and a = 0.5 in the case of 
the stochastic thermostat. The dashed lines are the Maxwell- 
Boltzmann predictions. 

4.2 

Gaussian thermostat 

Now we carry out a parallel analysis in the case of the 
deterministic Gaussian thermostat. The a-dependence of 
the simulation values of a-2, (n2 ~ A^Vm^ anc ^ (M4 — 
At4 '')/M4 are shown in Figure ^. The values of a 2 are 
in this case generally larger than in the previous case. In 
addition, Eq. (|2^) tends to overestimate /X4 and underes- 
timate fj,2 for small a. As a consequence, the theoretical 
estimate (JsJ) gives values larger than the simulation data 
for a < 0.5, while the estimate ( |40| ) is fairly good in that 
region. For values of the coefficient of restitution for which 
the fourth cumulant 02 is not small enough (say 02 — 0.1), 
we may expect a non-negligible deviation from (p4j). This 
is confirmed in Fig. [|, where A(c) is plotted for a = 0.4. 
Here the contributions associated with higher-order So- 
nine polynomials are relatively important. As a quantita- 
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0.20 - 



0.15 - 



0.10 - 



0.05 



0.00 - 



-0.05 



Fig. 4. Plot of the simulation values of a2 (O)j (a*2— l^)/!^ 
(A) and (/m — /U4 )//4, (v) versus a in the case of the Gaus- 
sian thermostat. The solid and dashed lines are the theoretical 
estimates (|B|) and (po|), respectively. 
30 




Fig. 5. Plot of the simulation values of the function A(c) de- 
fined by Eq. (^) for a = 0.4 in the case of the Gaussian 
thermostat. The dashed line is the Sonine polynomial 5*2(c 2 ). 

tive measure of the difference between A{c) and S2 (c 2 ), 
we have obtained preliminary simulation results for the 
sixth cumulant 03 defined as 
48 



03 



d{d + 2){d + A) 



d(d + 2)(d + 4) 



dcS 3 (c 2 )f(c) 



(c 6 ) + 1 + 3a 2 . 



Fig. 6. Plot of the simulation values of 0,3 versus a in the case 
of the Gaussian thermostat. 

This quantity is plotted in Fig. |^ for d — 3. For a > 0.6 
1 03) remains small, but for larger dissipation the values of 
I a.3 1 increase rapidly. 

The high energy tail predicted by Eq. ( |42] ) ||D is 
tested in Fig. where 



G(c) ^ e Ac f{c) 



(71) 



(70) 



is plotted for a = 0.2 and a = 0.4. The corresponding 
values of A are A ~ 3.82 and A ~ 4.41, respectively. The 
agreement with Eq. (|6^) is excellent; from the simulation 
data we can estimate K ~ 7 for a = 0.2 and K ~ 31 
for a = 0.4. In this case of a Gaussian thermostat, the 
overpopulation effect is much more important than in the 
previous case. At c = 4, f /<f) ~ 80 for a = 0.2 and f /<f) ~ 
34 for a = 0.4. The results reported here for inelastic 
hard spheres complement those obtained by Brey et al. 
[||, where the asymptotic behavior (^) was verified for 
inelastic hard disks. 

Recently, Sela and Goldhirsch |2(J have obtained nu- 
merically the function A{c) in the low dissipation limit. In 
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Fig. 7. Plot of the simulation values of the function G(c) de- 
fined by Eq. (fn] ) for a = 0.2 and a = 0.4 in the case of 
the Gaussian thermostat. The dashed lines are the Maxwell- 
Boltzmann predictions. 

their notation, lim Q ^i A(c) = — 8<? e (c). From simulation 
results presented in Ref. || for a = 0.99 it follows that 
the function <? e (c) is well represented by the Sonine poly- 
nomial — S , 2(c 2 )/8 in the range < c < 1. However, this 
agrees only qualitatively with the function <P C (c) obtained 
numerically by Sela and Goldhirsch | p0| . For instance, 
from Fig. 3 of Ref. |(J one gets $ e (0) ~ -0.35, while 
-S 2 (0)/8 = -15/64 ~ -0.23. Moreover, it is claimed in 
Ref. |2(J that <P e (c) ~ c 2 log c for large c, which differs 
from the behavior ( [i"2"| ) that has been confirmed here and 
in Ref. |j] . It is possible that the high energy tail obtained 
from the perturbative approach presented in Ref. [ p0[ only 
holds for 1 <C c <C (1 — a 2 ) -1 and thus it is not represen- 
tative of the general asymptotic behavior for arbitrary a. 



Fig. 8. Plot of the simulation values of 02 (0)> (/te - a4°') /Ma 
(A) and (/14 — /i^ ')/^^ 1 ' (v) versus a in the case of the non- 
Gaussian thermostat. The solid and dashed lines are the esti- 
mates ( |50| ) and (|T6j) , respectively. 

4.3 

Non-Gaussian thermostat 

In contrast to the two previous cases, the Sonine polyno- 
mials {S p (c 2 )} are not expected to constitute a good set 
for the expansion of the ratio /(c) /0(c) in the case of the 



non-Gaussian thermostat (15) since the latter is singular 
at c = 0. Consequently, we do not expect the estimate (|5C 
to be quantitatively accurate. This is confirmed in Fig. 



where we observe that Eq. (50) gives values that are about 
20% smaller in magnitude than the simulation ones. Also, 
the approximation ( |26| ) with ^\ given by Eqs. ( p8| ) and 
( |30| ) is rather poor. It is reasonable to expect that a better 
approximation would be obtained if (c), (c 3 ), [ii and /X4 
were computed from the unknown function A(c) rather 
than from S^c 2 ). When plotting the simulation data of 

(ya-yfip)/fi2> (M4-M4 \ ( c ) and ( c3 ) versus a 2 we 

have observed that the points fit well in straight lines, as 
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predicted by Eqs. (|26|), (|8|) and (|l9j), but with different 
slopes. More specifically, our simulation results indicate 
that, instead of Eqs. (p6|), ( [48| ) and (|49|), one should have 
(for d = 3) 

21 



(o) . ^ (1) 
- M 2 + ^2 °2, 



~ (0) , 13 (1) 

15 



2a/2 / 23 3 
^l 1 + 20 8 a2 



(72) 



(73) 



(74) 



(75) 



If we insert the above expressions into Eq. (44 ) and neglect 
terms nonlinear in 02, we get 

240(1 -a)(l + 2a 2 ) 



a.2 



(76) 



1957 - 1125a + 390(1 - a)a 2 ' 

This semi-empirical estimate exhibits a fairly good agree- 
ment with the simulation data, as shown in Fig. |8[ 

The limitations of a Sonine description in the case of 
the non-Gaussian thermostat are quite apparent in Fig. ^], 
where A(c) is plotted for a = 0.4, 0.6 and 0.95. The curves 
corresponding to a — 0.4 and a — 0.6 practically coincide, 
while the curve corresponding toa = 0.95 clearly deviates 
in the region of very small velocities. As a matter of fact, 
Z\(0) is roughly equal to — a? 1 , which indicates an almost 
vanishing population of rest particles, i.e. /(0) ~ 0, even 
at a = 0.95. A key feature of Fig. |9| is the existence of a 
non-zero initial slope, dA(c)/dc\ c=0 ^ 0, that cannot be 
described by any polynomial in c 2 . 



From the analysis made at the end of Subsect. 2.3 



we expect an underpopulated high energy tail of the form 
([S4|), where the coefficient A is unknown. By a fitting of 




Fig. 9. Plot of the simulation values of the function A(c) de- 
fined by Eq. @ for a = 0.4, 0.6 and 0.95 in the case of the 
non-Gaussian thermostat. The dashed line is the Sonine poly- 
nomial 52 (c 2 ). 

the simulation results we have estimated A ~ 1.48 for 
a = 0.3 and A ~ 1.51 for a — 0.4. Figure [l(] shows the 
function 



G(c) = e Ac -f(c) 



(77) 



for a = 0.3 and a = 0.4. In both cases the value of 
K is K ~ 1.7. The regions of small and large velocities 
are highly underpopulated with respect to the Maxwell- 
Boltzmann distribution. At c = 3, for instance, f /(f) — 0.12 
for a = 0.3 and f/(f> ~ 0.10 for a = 0.4. 



Summary and discussion 

In this paper we have performed direct Monte Carlo sim- 
ulations of the Enskog-Boltzmann equation for a fluid of 
smooth inelastic spheres in spatially uniform states. Upon 
describing the velocity distribution of the granular fluid 
by the Enskog-Boltzmann equation (|l|) it has been im- 



1 — i — 1 — i — 1 — i — 1 — i — 1 — i — 1 — i — : 

cc=oa ' X 




0.0 0.5 1.0 1.5 2.0 2.5 3.0 



c 

Fig. 10. Plot of the simulation values of the function G(c) 
defined by Eq. ( fTi] ) for a = 0.3 and a = 0.4 in the case of the 
non-Gaussian thermostat. The dashed lines are the Maxwell- 
Boltzmann predictions. 

plicitly assumed the validity of the "molecular chaos" hy- 
pothesis of uncorrelated binary collisions. However, molec- 
ular dynamics simulations of hard disks have shown a 
non-uniform distribution of impact parameters for high 
enough dissipation (a < 0.8) [|lj. In addition, there exist 
long range spatial correlations in density and flow fields 
which cannot be understood on the basis of the Enskog- 
Boltzmann equation [B3] . These two effects are associated 
with the appearance of the so-called cluster instability |^3| 
for systems sufficiently large. Since we have simulated di- 
rectly the spatially uniform equation (0), such an insta- 
bility is precluded in the simulations. 

To compensate for cooling effects associated with the 
inelasticity of collisions, three types of "thermostatting" 
external driving forces have been considered. We have an- 
alyzed the deviation of the steady-state velocity distribu- 
tion function from the Maxwell-Boltzmann distribution, 
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as measured by the fourth cumulant a-i and by the high 
energy tail. 

A simple mechanism for thermostatting the system 
is to assume that the particles are subjected to random 
kicks fl2|| , what mimics the effects of shaking or vibrat- 
ing the vessel [Q. If this stochastic force has the prop- 
erties of a white noise [cf. Eq. (||)], it gives rise to a 
Fokker-Planck diffusion term in the Enskog-Boltzmann 
equation ||. By making a first Sonine approximation, 
van Noije and Ernst g have obtained an approximate 
expression for ct2 as a function of the coefficient of normal 
restitution a. Our simulation results confirm the accuracy 
of that expression even for large dissipation (a = 0.2). 
We have also confirmed a high energy tail of the form 
/(v) ~ exp[— A(v/vo) 3 / 2 ] (where vq is the thermal ve- 
locity) derived in Ref. g. Moreover, that asymptotic be- 
havior (which represents an overpopulation with respect 
to the Maxwell-Boltzmann distribution) is already practi- 
cally reached for v > 4wo, at least for a = 0.4 and 0.5. 

In the absence of any external forcing, the freely evolv- 
ing granular fluid reaches a homogeneous cooling state in 
which all the time dependence of the velocity distribution 
occurs through the thermal velocity vq (t) , so that the dis- 
tribution of the reduced velocity c = v/vq is stationary. 
When the Enskog-Boltzmann equation is written in terms 
of this reduced velocity, the operator d/dt gives rise to 
an operator that coincides with the one representing the 
action of an external force proportional to the particle ve- 
locity [cf. Eq. (p^)]. This type of "anti-drag" force can also 
be justified by Gauss's principle of least constraint Jll| and 
has been widely used in nonequilibrium molecular dynam- 
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ics simulations of molecular fluids. Thus, the homogeneous and has now been confirmed by our simulation results for 
cooling state is equivalent to the steady state reached un- d = 3. 



der a Gaussian thermostat. In their simulations, Brey et 
al. |^,^] used the former point of view, while in this pa- 
per we have used the latter. Our simulations complement 
those of Ref. Q also in that we have considered a wide 
range 0.2 < a < 1, while Brey et al. (f| analyzed in de- 
tail the region 0.7 < a < 1. They obtained an excellent 
agreement with the estimate (|^) based on a Sonine approx- 
imation, first derived in Ref. Q|. However, as a decreases 
and a 2 becomes larger, we have seen in this paper that 
Eq. (j^) overestimates a-i- This discrepancy can be traced 
back to contributions associated with higher order Sonine 
polynomials as well as to the ambiguity involved in the 
approximate determination of ai by neglecting nonlinear 
terms in the exact equation (p8|). If Eq. ( |38| ) is rewritten 
in another equivalent form (for example, by transferring 
a quantity from one side to the other), the same method 
yields a different approximation for a-i- As long as a 2 re- 
mains small (say \ai\ < 0.05), all the approximations give 
practically undistinguishable results. On the other hand, 
for larger values of a-i (i.e., for a < 0.5) the result is rel- 
atively dependent of the route followed. By starting from 
Eq. ( |38| ) rewritten as 1^4/(1 + 02) = (<i + 2)/i2, we have ob- 



tained the estimate (40), which is seen to agree fairly well 
with the simulation results for the whole range of coeffi- 
cients of restitution considered. The asymptotic analysis of 
the kinetic equation predicts a high energy tail of the form 
HD /( v ) ~ ex P[~ ^W^o)], what represents an overpop- 
ulation phenomenon stronger than in the previous case. 
This behavior was already confirmed in Ref. (|] for d = 2 



In the case of the Gaussian thermostat, the heating 
force points in the motion direction and its magnitude is 
proportional to that of the particle velocity. This is a very 
efficient thermostat because it gives more energy to fast 
particles, which are the ones colliding more frequently. In 
contrast, the stochastic thermostat adds a velocity incre- 
ment per unit of time that is random both in direction 
and in magnitude. This is why the high energy popula- 
tion is larger with the Gaussian thermostat than with the 
stochastic thermostat. Nevertheless, in both cases such a 
population is larger than in the case of elastic particles 
at equilibrium. One could be tempted to expect that this 
overpopulation is a common feature of heated granular 
fluids, regardless of the mechanism of heating. Our third 
choice of thermostat, Eq. (jlj), proves that this is not the 
case. Like in the case of the stochastic thermostat, the 
force is independent of the magnitude of the particle ve- 
locity; like in the case of the Gaussian thermostat, the 
force is deterministic and points in the motion direction. 
The action of this third thermostat can be graphically de- 
scribed by saying that, between two successive collisions, 
a particle feels a "pseudo-gravity" field that makes it to 
"fall" along its motion direction. With this choice of a 
non-Gaussian deterministic thermostat, the Sonine poly- 
nomials {S p (c 2 )} are not a good set to represent the ra- 
tio /(c)/0(c), even for low dissipation. As a consequence, 
the theoretical estimate of 02 derived by assuming that 
[/(c) /0(c) - l]/o 2 = A{c) ~ 5 2 (c 2 ), while being quali- 
tatively correct, is not quantitatively accurate. We have 
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not been able to get the functional form of A(c) in the 
limit of low dissipation. However, we have estimated its 
contributions to (c), (c 3 ), [ii and /14 from the simulation 
data. This has allowed us to obtain an approximate ex- 
pression for d2 that fits well the simulation results. An 
interesting feature of the velocity distribution function in 
this case is that it is highly underpopulated with respect 
to the Maxwell-Boltzmann distribution both for small and 
large velocities. Between two successive collisions, every 
particle experiences a constant tangential acceleration g. 
The total work done by this force is exactly compensated 
by the total loss of energy through collisions, which are 
much more frequent for fast particles than for slow ones. 
Therefore, the population of slow particles decreases be- 
cause of the action of the external force, while that of fast 
particles decreases because of the effect of collisions. The 
high energy tail of the distribution function is of the form 
/(c) ~ exp(— Ac 2 ) with A > 1. In this case the gain and 
loss terms of the collision integral are comparable, so that 
the dependence of A on a is an open problem. 

Partial support from the DGES (Spain) through grant No. 
PB97-1501 and from the Junta de Extremadura-Fondo Social 
Europeo through grant No. IPR98C019 is gratefully acknowl- 
edged. 



it is easy to get 



Collisional moments 



fi p = J dcx J dc 2 /(ci)/(c2)#p(ci,C2), (78) 



where 



# p (ci,C2) = - / da 0(c 12 ■ ff)(cia ■ a) [c? + c§ - c'/ - c' 2 p ] 

(79) 

with c[ 2 = ci j2 =F |(1 + cc)(ci2 ■ a)a. In the cases p — 2 
and p — 4 we have 

I -a 2 



2,2 ' 2 ,2 
c l T c 2 — c l — c 2 



-{c l2 -a) 2 , (80) 



'■i+r.j-r, -r L , = (1 - a 2 ) (c 12 • a) 2 [ -±* + C 



(1-a 2 ) 2 



g (ci2-^) 4 

-2(l + a) 2 (c 12 -S) 2 (C 12 -£) 2 
+4(l + a)(ci2 -ct)(Ci2 -a) 
x(c 12 -C 12 ), (81) 
where Ci 2 = c± — c 2 , C12 = §(ci + C2). Consequently, 

(82) 



5C / n /3 3 (l-a 2 ) 3 
^ 2 (ci,c 2 ) = £-i— J -c 3 



12- 



^4(C1,C 2 ) = /3 3 Ci2 



(d + 2)(l-a 2 ) a 

2d 12 12 



2 ,-<2 



(l-q 2 )(rf+l + 2q 2 ) 4 
8(d + 3) ° 12 
(2d + 3-3a)(l +a) 
d + 3 

(C12 • Ci2) 2 — ^C 2 2 C 2 2 



where we have taken into account that 



(83) 



J da 0(c ■ a)(c ■ a) n a — (3 n +ic 



da 6>(c ■ a)(c ■ a) n aa 



n + d 

In this Appendix we derive the expressions (|6^)-(|66|). Start- In the three-dimensional case, Eqs. (| 
ing from Eq. ( Egh and by a standard change of variables, yield Eqs. (^)-(^6|). 



(ncc + I) 



(84) 



(85) 



(86) 



3) and (| 
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